#include "DynamicLoad.h"

#include "DataStructure/Common/Constant.h"
#include "DataStructure/Input/Solver/Subcase.h"
#include "DataStructure/Input/Table/TABLED1Data.h"


#include "DBManager/IService/Input/Load/IDAREAService.h"
#include "DBManager/IService/Input/Load/IDELAYService.h"
#include "DBManager/IService/Input/Subcase/ISubcaseService.h"
#include "DBManager/IService/Input/Table/ITABLED1Service.h"

#include "SGFem/Pre/Assembly/Assembly.h"

#include "SGFem/Task/Common/LoadCalculator.h"

using SG::Algebra::Real_t;
using SG::Algebra::Complex_t;
using namespace SG::DataStructure::Common;
using namespace SG::DataStructure::FEM;
using namespace SG::DBManager;
using SG::FEMSolver::AssemblyEigen;

/**
 * @brief    计算空间载荷分布
 *
 * @param    db
 * @param    subcaseId 工况号，静力载荷计算查询功能
 * @param    loadsetId loadset 具有loadset 时查询 lseq
 * @param    excitId   当具有lseq 时，和 loadset 一起查询 lseq; 否则用于查询静力载荷
 * @param[out]    fVec 约束处理完毕的载荷右端项
 * @author   wanghuan
 * @date     2025.09.08
 */
void calcLoad (SG::DBManager::DBServiceFactorySP db,
               SG::DataStructure::Common::Id     subcaseId,
               Id                                loadsetId,
               SG::DataStructure::Common::Id     excitId,
               SG::Algebra::Vec_t&               fVec)
{
    auto& assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    // DAREA
    std::vector<SG::DataStructure::FEM::DAREAData> DAREAList;
    if (Status::STATUS_Success == db->get<IDAREAService> ()->FindBySId (excitId, DAREAList))
    {
        LoadCalculator::calcDAREAForce (DAREAList, fVec);
    }

    SG::Algebra::Vec_t fStatic;

    // Static Load
    Id loadId = excitId;
    // has LoadSet ?
    if (0 != loadsetId)
    {
        // loadId = LSEQ.m_LID
        SG::DataStructure::FEM::LSEQData lseq;
        db->get<SG::DBManager::ILSEQService> ()->FindByLoadId (loadsetId, excitId, lseq);
        loadId = lseq.m_LID;
    }

    // get load vector
    LoadCombination load;
    const auto      isLoadCase = db->get<ILoadCombinationService> ()->FindById (loadId, load);
    if (Status::STATUS_Success == isLoadCase)
    {
        assembly.ComputeBVector (db, load, _OUT fStatic);
    }
    else
    {
        assembly.ComputeBVector (loadId, _OUT fStatic);
    }

    fVec += fStatic;
}

SG::DataStructure::FEM::TABLEDObjectP registTimeTable (const SG::DBManager::DBServiceFactorySP db, const Id tId)
{
    Status status = Status::STATUS_Failure;
    
    TABLED1Data tabled1;
    status = db->get<ITABLED1Service> ()->FindById (tId, tabled1);
    if (Status::STATUS_Success == status)
    {
        return std::make_shared<TABLED1Data>(tabled1);
    }

    return NULL;
}


TLoad1Calc::TLoad1Calc (const SG::DataStructure::FEM::TLOAD1Data& tload1,
                        SG::DataStructure::Common::Id             loadsetId,
                        SG::DataStructure::Common::Real           S)
    : m_load (tload1), m_timeTab {}, m_loadsetId (loadsetId), m_S (S)
{}

SG::DataStructure::Common::Status TLoad1Calc::computeA (SG::DBManager::DBServiceFactorySP db, const SG::DataStructure::Common::Id subcaseId)
{
    auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    std::size_t n        = assembly.GetFreeDOFSize ();
    m_staticF            = SG::Algebra::Vec_t (n);
    m_staticF.setZero ();

    const auto excitId = m_load.m_DAREA;
    calcLoad (db, subcaseId, m_loadsetId, excitId, _OUT m_staticF);

    if (m_load.m_TID > 0)
    {
        m_timeTab = registTimeTable (db, m_load.m_TID);
    }

    // DELAY
    m_delayMap.clear ();
    if (m_load.m_DELAY > 0)
    {
        std::vector<SG::DataStructure::FEM::DELAYData> delayList;
        if (Status::STATUS_Success == db->get<IDELAYService> ()->FindBySId (m_load.m_DELAY, _OUT delayList))
        {
            for (const auto& delay : delayList)
            {
                std::size_t nodeIndex;
                int         dofGlobalId;
                assembly.FindNodeIndex (delay.m_P, _OUT nodeIndex);
                assembly.GetGlobalDOFId (nodeIndex, delay.GetDirection (), _OUT dofGlobalId);
                m_delayMap[dofGlobalId] = delay.m_T;
            }
        }
    }

    return SG::DataStructure::Common::Status::STATUS_Success;
}

SG::Algebra::Vec_t TLoad1Calc::Interplate (SG::DBManager::DBServiceFactorySP db, const SG::Algebra::Real_t t) const
{
    SG::Algebra::Vec_t result (m_staticF.size ());
    result.setZero ();

    // 时间插值
    auto globalCoeff = m_load.m_F;
    if (m_load.m_TID > 0)
    {
        // 所有自由度使用一致的延时DELAYR，默认0.0
        const auto tau = m_load.m_T;
        globalCoeff = m_timeTab -> Interpolate (t - tau);
    }

    result           = (globalCoeff * m_S) * m_staticF;

    // 载荷随时间变化（即用插值表），且某自由度具有指定的延时（使用DELAYI），需要重新计算对应自由度载荷值
    if (m_load.m_TID > 0 && m_load.m_DELAY > 0)
    {
        // 根据DELAY词条对存在单独延迟的自由度做计算覆盖
        for (const auto& [dof, delay] : m_delayMap)
        {
            SG::Algebra::Real_t coeff = m_timeTab -> Interpolate (t - delay);
            result[dof]               = coeff * m_S * m_staticF[dof];
        }
    }

    return result;
}

TLoad2Calc::TLoad2Calc (const SG::DataStructure::FEM::TLOAD2Data& tload2,
                        SG::DataStructure::Common::Id             loadSetId,
                        SG::DataStructure::Common::Real           S)
    : m_load (tload2), m_loadsetId (loadSetId), m_S (S)
{}

SG::DataStructure::Common::Status TLoad2Calc::computeA (SG::DBManager::DBServiceFactorySP db, const SG::DataStructure::Common::Id subcaseId)
{
    auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    std::size_t n        = assembly.GetFreeDOFSize ();
    m_staticF            = SG::Algebra::Vec_t (n);
    m_staticF.setZero ();
    const auto excitId = m_load.m_DAREA;
    calcLoad (db, subcaseId, m_loadsetId, excitId, _OUT m_staticF);

    // DELAY
    m_delayMap.clear ();
    if (m_load.m_DELAY > 0)
    {
        std::vector<SG::DataStructure::FEM::DELAYData> delayList;
        if (Status::STATUS_Success == db->get<IDELAYService> ()->FindBySId (m_load.m_DELAY, _OUT delayList))
        {
            for (const auto& delay : delayList)
            {
                std::size_t nodeIndex;
                int         dofGlobalId;
                assembly.FindNodeIndex (delay.m_P, _OUT nodeIndex);
                assembly.GetGlobalDOFId (nodeIndex, delay.GetDirection (), _OUT dofGlobalId);
                m_delayMap[dofGlobalId] = delay.m_T;
            }
        }
    }

    return SG::DataStructure::Common::Status::STATUS_Success;
}

SG::Algebra::Vec_t TLoad2Calc::Interplate (SG::DBManager::DBServiceFactorySP db, const SG::Algebra::Real_t t) const
{
    SG::Algebra::Vec_t result (m_staticF.size ());
    result.setZero ();

    auto Calccoff = []  (const TLOAD2Data& load,Real_t tHat) -> Real_t
    {
        return std::pow (tHat, load.m_B) * std::exp (load.m_C * tHat) * cos (2 * PI * (load.m_F * tHat + load.m_P / 360));;
    };

    // DELAY是实数 && 没有DELAY
    SG::Algebra::Real_t uniformCoeff = 0.0;
    if (t >= (m_load.m_T1 + m_load.m_T) && t <= (m_load.m_T2 + m_load.m_T))
    {
        auto tHat = t - m_load.m_T1 - m_load.m_T;
        uniformCoeff = Calccoff (m_load, tHat);     
    }
    result = (uniformCoeff * m_S) * m_staticF;

    // 根据DELAY词条对存在单独延迟的自由度做计算覆盖
    for (const auto& [dof, delay] : m_delayMap)
    {
        SG::Algebra::Real_t coeff = 0.0;
        if (t >= (m_load.m_T1 + delay) && t <= (m_load.m_T2 + delay))
        {
            auto tHat = t - m_load.m_T1 - delay;
            coeff = Calccoff (m_load, tHat);
        }
        result[dof] = coeff * m_S * m_staticF[dof];
    }

    return result;
}

RLoad1Calc::RLoad1Calc (const SG::DataStructure::FEM::RLOAD1Data& rload1, Id loadsetId, SG::DataStructure::Common::Real S)
    : m_rload1 (rload1), m_TC (), m_TD (), m_S (S), m_loadsetId (loadsetId)
{}

SG::DataStructure::Common::Status RLoad1Calc::computeA (SG::DBManager::DBServiceFactorySP db, const SG::DataStructure::Common::Id subcaseId)
{
    auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    const std::size_t n        = assembly.GetFreeDOFSize ();
    m_staticF            = SG::Algebra::Vec_t (n);
    m_staticF.setZero ();

    if (DynaExcitType::DYNAEXC_Force == m_rload1.GetExcitType ())
    {
        const auto excitId = m_rload1.m_DAREA;
        calcLoad (db, subcaseId, m_loadsetId, excitId, _OUT m_staticF);
    }
    
    // 查询表
    if (m_rload1.m_TC > 0)
    {
        m_TC = registTimeTable (db, m_rload1.m_TC);
    }

    if  (m_rload1.m_TD > 0)
    {
        m_TD = registTimeTable (db, m_rload1.m_TD);
    }

    //延迟
    if (m_rload1.m_DELAY > 0)
    {

    }

    //相位
    if (m_rload1.m_DPHASE > 0)
    {
        
    }

    return SG::DataStructure::Common::Status::STATUS_Success;
}

SG::Algebra::ComplexVec_t RLoad1Calc::Interplate (SG::DBManager::DBServiceFactorySP db, const SG::Algebra::Real_t t) const
{
    auto C = m_rload1.m_RC;
    if (m_rload1.m_TC > 0 && NULL != m_TC)
    {
        C = m_TC -> Interpolate (t);
    }

    auto D = m_rload1.m_RD;
    if (m_rload1.m_TD > 0 && NULL != m_TD)
    {
        D = m_TD -> Interpolate (t);
    }
    auto coeff = Complex_t (C, D);

    // ToDo 全局的相位和延迟处理
    // Real_t theta = m_rload1.m_PH - 2*PI*t*m_rload1.m_T;
    const auto exitType = m_rload1.GetExcitType ();
    if (DynaExcitType::DYNAEXC_Force == exitType)
    {
        return (coeff * m_S) * m_staticF;
    }
    else
    {
        auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
        // spcd 被 load引用 不会发生？
        SG::Algebra::ComplexVec_t force;
        assembly.ComputeFreqSPCDLoad (m_rload1.m_DAREA, exitType, t, force);
        return (coeff * m_S) * force;
    }

    // ToDo 某个自由度上有相位差&延迟
    if(m_rload1.m_DELAY > 0 || m_rload1.m_DPHASE)
    {
        // 这些自由度的数值需要更新

    }
}

RLoad2Calc::RLoad2Calc (const SG::DataStructure::FEM::RLOAD2Data& rload2, Id loadsetId, SG::DataStructure::Common::Real S)
    : m_rload2 (rload2), m_timeTabB (), m_S (S), m_loadsetId (loadsetId)
{}

SG::DataStructure::Common::Status RLoad2Calc::computeA (SG::DBManager::DBServiceFactorySP db, const SG::DataStructure::Common::Id subcaseId)
{
    auto&       assembly = SG::FEMSolver::AssemblyEigen::GetInstance ();
    std::size_t n        = assembly.GetFreeDOFSize ();
    m_staticF            = SG::Algebra::Vec_t (n);
    m_staticF.setZero ();
    const auto excitId = m_rload2.m_DAREA;
    calcLoad (db, subcaseId, m_loadsetId, excitId, _OUT m_staticF);

    if (m_rload2.m_TB > 0)
    {
        m_timeTabB = registTimeTable (db, m_rload2.m_TB);
    }


    m_delayMap.clear ();
    if (m_rload2.m_DELAY > 0)
    {
        std::vector<SG::DataStructure::FEM::DELAYData> delayList;
        if (Status::STATUS_Success == db->get<IDELAYService> ()->FindBySId (m_rload2.m_DELAY, delayList))
        {
            for (const auto& delay : delayList)
            {
                std::size_t nodeIndex;
                int         dofGlobalId;
                assembly.FindNodeIndex (delay.m_P, _OUT nodeIndex);
                assembly.GetGlobalDOFId (nodeIndex, delay.GetDirection (), _OUT dofGlobalId);
                m_delayMap[dofGlobalId] = delay.m_T;
            }
        }
    }

    return SG::DataStructure::Common::Status::STATUS_Success;
}

SG::Algebra::ComplexVec_t RLoad2Calc::Interplate (SG::DBManager::DBServiceFactorySP db, const SG::Algebra::Real_t t) const
{
    SG::Algebra::ComplexVec_t result (m_staticF.size ());
    result.setZero ();

    double uniform_delay = m_rload2.m_T;

    auto global_B_f{m_rload2.m_RB};
    if (m_rload2.m_TB > 0)
    {
        global_B_f = m_timeTabB->Interpolate(t - uniform_delay);
    }

    result = (global_B_f * m_S) * m_staticF;

    for (const auto& [dof, delay] : m_delayMap)
    {
        double dof_B_f{m_rload2.m_RB};
        if (m_rload2.m_TB > 0)
        {
            dof_B_f = m_timeTabB->Interpolate(t - delay);
        }

        result[dof] = dof_B_f * m_S * m_staticF[dof];
    }
    
    return result;
}